library(here)library(dplyr)library(ggplot2)library(sf)library(rnaturalearth)library(smoothr)library(raster)library(readr)library(tidyr)library(ggpubr)library(betareg)library(broom)library(margins)library(patchwork)library(dplyr)library(sf)library(egg)library(spatstat)library(sf)library(tidyverse) # For data manipulation and visualizationlibrary(rnaturalearthdata) # Additional natural earth datalibrary(biscale) # For bivariate mappinglibrary(gridExtra) # For arranging multiple plotslibrary(grid) # For text elements in plotslibrary(colorspace) # For color manipulationlibrary(knitr)library(kableExtra)
Load the data
Code
# Read the CAPTAIN2 EDGE2 RDS fileCAPTAIN2_EDGE2_data <-readRDS(here::here("Data/CAPTAIN2_EDGE_full_results_averaged_budget0.1_replicates50.rds"))# Read the CAPTAIN2 FUSE RDS fileCAPTAIN2_FUSE_data <-readRDS(here::here("Data/CAPTAIN2_FUSE_res_full_results_averaged_budget0.1_replicates50.rds"))# Read the CAPTAIN2 IUCN RDS fileCAPTAIN2_IUCN_data <-readRDS(here::here("Data/CAPTAIN2_IUCN_full_results_averaged_budget0.1_replicates50.rds"))# Load one of your input raster files to extract the correct grid structure# Use the same raster file for consistencyraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Read the protected range fractions RDS fileprotected_fractions <-readRDS(here::here("Data", "CAPTAIN2_protected_range_fractions_2ndrun.rds"))# Read the continental shark conservation metrics CSV fileshark_metrics <-read_csv(here::here("Data", "continental_shark_conservation_metrics_10_harmonised_IUCN_categories.csv"))# Load fishing dataload(here::here("Data", "Raw", "Predicted_Fishing_Hours_05Deg.Rdata"))
Summary Statistics of IUCN Conservation Priority Values
Cell Counts
All Cells Statistics
Non-zero Only
Index
Total Cells
Non-zero Cells
Non-zero %
Zero Cells
Zero %
NA Cells
NA %
Min
Q25
Median
Mean
Q75
Max
Non-zero Median
Non-zero Mean
IUCN
232560
6031
2.59
226529
97.41
0
0
0
0
0
0.02
0
1
1
0.85
Code
# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 IUCN data based on PUIDCAPTAIN2_IUCN_data_with_coords <- CAPTAIN2_IUCN_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_IUCN_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 IUCN data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_IUCN_data_nonzero <- CAPTAIN2_IUCN_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_IUCN_sf <-st_as_sf( CAPTAIN2_IUCN_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_IUCN <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_IUCN <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_IUCN <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank(),axis.text =element_blank(), # Add this lineaxis.ticks =element_blank() # Add this line )# Create the plotCAPTAIN2_IUCN_plot <-ggplot() +geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_IUCN, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_IUCN, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_IUCN# Display the plotprint(CAPTAIN2_IUCN_plot)
Code
# Save the plotggsave(filename = here::here("Outputs", "CAPTAIN2_IUCN_priorities_01_2.png"),plot = CAPTAIN2_IUCN_plot,width =10,height =6,dpi =300,bg ="white")
Summary Statistics of FUSE Conservation Priority Values
Cell Counts
All Cells Statistics
Non-zero Only
Index
Total Cells
Non-zero Cells
Non-zero %
Zero Cells
Zero %
NA Cells
NA %
Min
Q25
Median
Mean
Q75
Max
Non-zero Median
Non-zero Mean
FUSE
232560
7271
3.13
225289
96.87
0
0
0
0
0
0.02
0
1
1
0.75
Code
# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cell# This gives us the center coordinates of each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 FUSE data based on PUIDCAPTAIN2_FUSE_data_with_coords <- CAPTAIN2_FUSE_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_FUSE_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 FUSE data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_FUSE_data_nonzero <- CAPTAIN2_FUSE_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_FUSE_sf <-st_as_sf( CAPTAIN2_FUSE_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_FUSE <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_FUSE <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_FUSE <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank(),axis.text =element_blank(), # Add this lineaxis.ticks =element_blank()) # Add this line# Create the plotCAPTAIN2_FUSE_plot <-ggplot() +geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_FUSE, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_FUSE, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_FUSE# Display the plotprint(CAPTAIN2_FUSE_plot)
Code
# Save the plotggsave(filename = here::here("Outputs", "CAPTAIN2_FUSE_priorities_01_2.png"),plot = CAPTAIN2_FUSE_plot,width =10,height =6,dpi =300,bg ="white")
Summary Statistics of EDGE2 Conservation Priority Values
Cell Counts
All Cells Statistics
Non-zero Only
Index
Total Cells
Non-zero Cells
Non-zero %
Zero Cells
Zero %
NA Cells
NA %
Min
Q25
Median
Mean
Q75
Max
Non-zero Median
Non-zero Mean
EDGE2
232560
5618
2.42
226942
97.58
0
0
0
0
0
0.02
0
1
1
0.96
Code
# Load one of your input raster files to extract the correct grid structureraster_file <- here::here("Data", "tif files continental", "Psammobatis_parvacauda.tif")# Check if the file existsif (!file.exists(raster_file)) {stop("Raster file not found. Please provide a valid path to one of your input raster files.")}# Load the rasterr <-raster(raster_file)# Get the dimensions of the rasternrows <-nrow(r)ncols <-ncol(r)# Confirm dimensions match expected valuesif (nrows !=323|| ncols !=720) {warning("Raster dimensions don't match expected values. Proceeding with actual dimensions.")}# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Now join with the CAPTAIN2 data based on PUIDCAPTAIN2_EDGE2_data_with_coords <- CAPTAIN2_EDGE2_data %>%left_join(coords, by ="PUID")# Check if the join worked correctlyif (sum(is.na(CAPTAIN2_EDGE2_data_with_coords$Longitude)) >0) {warning("Some PUIDs from CAPTAIN2 data couldn't be matched to coordinates.")}# Filter to keep only cells with non-zero priority for faster plottingCAPTAIN2_EDGE2_data_nonzero <- CAPTAIN2_EDGE2_data_with_coords %>%filter(Priority >0) %>%filter(!is.na(Longitude), !is.na(Latitude)) # Remove any rows with missing coords# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Transform the dataset to sf object and projectCAPTAIN2_EDGE2_sf <-st_as_sf( CAPTAIN2_EDGE2_data_nonzero, coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE) # Use the raster's CRS) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected_CAPTAIN2_EDGE2 <-st_transform(world, crs = mcbryde_thomas_2)# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border_CAPTAIN2_EDGE2 <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Create base thememy_theme_CAPTAIN2_EDGE2 <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank(),axis.text =element_blank(), # Add this lineaxis.ticks =element_blank() # Add this line )# Create the plotCAPTAIN2_EDGE2_plot <-ggplot() +geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Global Conservation Priorities",subtitle ="CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_EDGE2# Display the plotprint(CAPTAIN2_EDGE2_plot)
Code
# Save the plotggsave(filename = here::here("Outputs", "CAPTAIN2_2ndrun_EDGE2_priorities_01_2.png"),plot = CAPTAIN2_EDGE2_plot,width =10,height =6,dpi =300,bg ="white")
Difference maps
Code
# Create a grid of coordinates for each cellcoords <-as.data.frame(coordinates(r))names(coords) <-c("Longitude", "Latitude")# Add cell IDs (PUID) to the coordinatescoords$PUID <-1:nrow(coords)# Join all three datasets with coordinatesIUCN_with_coords <- CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority) %>%left_join(coords, by ="PUID")EDGE2_with_coords <- CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority) %>%left_join(coords, by ="PUID") FUSE_with_coords <- CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority) %>%left_join(coords, by ="PUID")# Combine all datasetsall_indices <- coords %>%left_join(CAPTAIN2_IUCN_data %>% dplyr::select(PUID, IUCN = Priority), by ="PUID") %>%left_join(CAPTAIN2_EDGE2_data %>% dplyr::select(PUID, EDGE2 = Priority), by ="PUID") %>%left_join(CAPTAIN2_FUSE_data %>% dplyr::select(PUID, FUSE = Priority), by ="PUID")# Calculate differencesall_indices <- all_indices %>%mutate(# Replace NA with 0 for calculation purposesIUCN =ifelse(is.na(IUCN), 0, IUCN),EDGE2 =ifelse(is.na(EDGE2), 0, EDGE2),FUSE =ifelse(is.na(FUSE), 0, FUSE),# Calculate differencesIUCN_minus_FUSE = IUCN - FUSE,IUCN_minus_EDGE2 = IUCN - EDGE2,EDGE2_minus_FUSE = EDGE2 - FUSE )# Get world map dataworld <-ne_countries(scale ="medium", returnclass ="sf")# Define the McBryde-Thomas 2 projectionmcbryde_thomas_2 <-"+proj=mbt_s"# Create the globe bounding boxglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))# Create the globe borderglobe_border <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Project the world mapworld_projected <-st_transform(world, crs = mcbryde_thomas_2)# Create base thememy_theme <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank(),axis.text =element_blank(), # Add this lineaxis.ticks =element_blank() # Add this line )# Filter to non-zero differences for each comparison to reduce plot size# IUCN - FUSEIUCN_FUSE_diff <- all_indices %>%filter(IUCN_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# IUCN - EDGE2IUCN_EDGE2_diff <- all_indices %>%filter(IUCN_minus_EDGE2 !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# EDGE2 - FUSEEDGE2_FUSE_diff <- all_indices %>%filter(EDGE2_minus_FUSE !=0) %>%st_as_sf(coords =c("Longitude", "Latitude"), crs =crs(r, asText =TRUE)) %>%st_transform(crs = mcbryde_thomas_2)# Create a diverging color palette for difference maps# Blue for negative (first index lower), white for zero, red for positive (first index higher)diff_colors <-c("#2166AC", "#4393C3", "#92C5DE", "#D1E5F0", "#FFFFFF", "#FDDBC7", "#F4A582", "#D6604D", "#B2182B")# 1. IUCN - FUSE Difference MapIUCN_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority (IUCN - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus FUSE Index",x =NULL, y =NULL) + my_theme# 2. IUCN - EDGE2 Difference MapIUCN_EDGE2_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority (IUCN - EDGE2)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="IUCN Index minus EDGE2 Index",x =NULL, y =NULL) + my_theme# 3. EDGE2 - FUSE Difference MapEDGE2_FUSE_plot <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="black", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority (EDGE2 - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(title ="Difference in Conservation Priorities",subtitle ="EDGE2 Index minus FUSE Index",x =NULL, y =NULL) + my_theme# Display all plotsprint(IUCN_FUSE_plot)
# Rename column in shark_metrics to match bettershark_metrics <- shark_metrics %>%rename(Species =`Species name`)# Order shark_metrics alphabetically by Species nameshark_metrics <- shark_metrics %>%arrange(Species)# Add an original order ID to protected_fractions to maintain its original orderprotected_fractions$original_order <-1:nrow(protected_fractions)# Add row number as IDs to both datasetsprotected_fractions$protected_ID <-1:nrow(protected_fractions)shark_metrics$species_ID <-1:nrow(shark_metrics)# Check if the datasets have the same number of rowsif(nrow(protected_fractions) ==nrow(shark_metrics)) {# Create index mapping indices <-data.frame(protected_ID =1:nrow(protected_fractions),species_ID =1:nrow(shark_metrics) )# Join protected_fractions with indices protected_with_indices <- protected_fractions %>%left_join(indices, by ="protected_ID")# Join shark_metrics with indices shark_with_indices <- shark_metrics %>%left_join(indices, by ="species_ID")# Now join the datasets combined_data <- protected_with_indices %>%inner_join( shark_with_indices,by =c("species_ID", "protected_ID"),suffix =c("_captain", "_original") ) %>%arrange(original_order)cat("Successfully joined datasets with", nrow(combined_data), "species\n")# Define IUCN categories and order iucn_order <-c("LC", "NT", "VU", "EN", "CR") iucn_colors <-c("LC"="#50C878","NT"="#FFFF00","VU"="#FFA500","EN"="#FF8C00","CR"="#FF0000" )# Calculate sample sizes for each group iucn_n <- combined_data %>%group_by(IUCN_original) %>%summarise(n =n()) %>%arrange(IUCN_original) fuse_n <- combined_data %>%group_by(FUSE_original) %>%summarise(n =n()) %>%arrange(FUSE_original) edge2_n <- combined_data %>%group_by(EDGE2_original) %>%summarise(n =n()) %>%arrange(EDGE2_original)# Create x-axis labels with sample sizes iucn_labels <-paste0(iucn_order, "\n(n=", iucn_n$n, ")") fuse_labels <-paste0(1:5, "\n(n=", fuse_n$n, ")") edge2_labels <-paste0(1:5, "\n(n=", edge2_n$n, ")")# 1. IUCN plot iucn_plot <- combined_data %>%mutate(IUCN_status =factor(IUCN_original, levels =1:5, labels = iucn_order),protection_percentage = IUCN_captain *100 ) %>%ggplot(aes(x = IUCN_status, y = protection_percentage)) +geom_jitter(width =0.1,size =0.4,alpha =0.5,color ="darkgray") +stat_summary(fun = mean,geom ="point",size =3,color ="black") +stat_summary(fun.data =function(x) { mean_val <-mean(x, na.rm =TRUE) sd_val <-sd(x, na.rm =TRUE)return(data.frame(y = mean_val,ymin =max(0, mean_val - sd_val),ymax =min(100, mean_val + sd_val))) },geom ="errorbar",width =0.1,color ="black") +labs(x ="IUCN Red List threat status",y ="Range protected (%)" ) +theme_minimal() +theme(plot.title =element_text(hjust =0.5, face ="bold", size =18),legend.position ="none",panel.grid.major.x =element_blank(),plot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="white", color =NA),axis.title.x =element_text(size =16),axis.title.y =element_text(size =16),axis.text.x =element_text(size =14),axis.text.y =element_text(size =14) ) +scale_x_discrete(limits = iucn_order, labels = iucn_labels) +scale_fill_manual(values = iucn_colors) +scale_color_manual(values = iucn_colors) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 2. FUSE plot fuse_plot <- combined_data %>%mutate(FUSE_category =factor(FUSE_original),protection_percentage = FUSE_captain *100 ) %>%ggplot(aes(x = FUSE_category, y = protection_percentage)) +geom_jitter(width =0.1,size =0.4,alpha =0.5,color ="darkgray") +stat_summary(fun = mean,geom ="point",size =3,color ="black") +stat_summary(fun.data =function(x) { mean_val <-mean(x, na.rm =TRUE) sd_val <-sd(x, na.rm =TRUE)return(data.frame(y = mean_val,ymin =max(0, mean_val - sd_val),ymax =min(100, mean_val + sd_val))) },geom ="errorbar",width =0.1,color ="black") +labs(x ="FUSE Score",y ="Range protected (%)" ) +theme_minimal() +theme(plot.title =element_text(hjust =0.5, face ="bold", size =18),legend.position ="none",panel.grid.major.x =element_blank(),plot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="white", color =NA),axis.title.x =element_text(size =16),axis.title.y =element_text(size =16),axis.text.x =element_text(size =14),axis.text.y =element_text(size =14) ) +scale_x_discrete(labels = fuse_labels) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# 3. EDGE2 plot edge2_plot <- combined_data %>%mutate(EDGE2_category =factor(EDGE2_original),protection_percentage = EDGE2_captain *100 ) %>%ggplot(aes(x = EDGE2_category, y = protection_percentage)) +geom_jitter(width =0.1,size =0.4,alpha =0.5,color ="darkgray") +stat_summary(fun = mean,geom ="point",size =3,color ="black") +stat_summary(fun.data =function(x) { mean_val <-mean(x, na.rm =TRUE) sd_val <-sd(x, na.rm =TRUE)return(data.frame(y = mean_val,ymin =max(0, mean_val - sd_val),ymax =min(100, mean_val + sd_val))) },geom ="errorbar",width =0.1,color ="black") +labs(x ="EDGE2 Score",y ="Range protected (%)" ) +theme_minimal() +theme(plot.title =element_text(hjust =0.5, face ="bold", size =18),legend.position ="none",panel.grid.major.x =element_blank(),plot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="white", color =NA),axis.title.x =element_text(size =16),axis.title.y =element_text(size =16),axis.text.x =element_text(size =14),axis.text.y =element_text(size =14) ) +scale_x_discrete(labels = edge2_labels) +scale_y_continuous(limits =c(0, 100),breaks =seq(0, 100, 25))# Combine the three plots into a grid combined_plots <- ggpubr::ggarrange( iucn_plot, fuse_plot, edge2_plot,ncol =3,nrow =1,labels =c("(A)", "(B)", "(C)"),font.label =list(size =14,face ="bold",hjust =-2,vjust =1.5) )#print(combined_plots)# Save combined plotggsave(here::here("Outputs", "Fig.2_combined_protection_dotplots.png"), combined_plots, width =20, height =6, dpi =150, bg ="white")#======================# Beta Regression Analysis#======================# Prepare data for beta regression model_data <- combined_data %>%mutate(IUCN_status = IUCN_original,FUSE_score = FUSE_original,EDGE2_score = EDGE2_original,IUCN_protection = IUCN_captain *100,FUSE_protection = FUSE_captain *100,EDGE2_protection = EDGE2_captain *100,# Convert to (0,1) scale and handle boundariesn_obs =n(),IUCN_protection_beta = (IUCN_protection/100* (n_obs -1) +0.5) / n_obs,FUSE_protection_beta = (FUSE_protection/100* (n_obs -1) +0.5) / n_obs,EDGE2_protection_beta = (EDGE2_protection/100* (n_obs -1) +0.5) / n_obs )# Fit beta regression models iucn_beta <-betareg(IUCN_protection_beta ~ IUCN_status, data = model_data) fuse_beta <-betareg(FUSE_protection_beta ~ FUSE_score, data = model_data) edge2_beta <-betareg(EDGE2_protection_beta ~ EDGE2_score, data = model_data)# View summariescat("\n=== IUCN Beta Regression ===\n")print(summary(iucn_beta))cat("\n=== FUSE Beta Regression ===\n")print(summary(fuse_beta))cat("\n=== EDGE2 Beta Regression ===\n")print(summary(edge2_beta))# Get tidy summaries iucn_tidy <-tidy(iucn_beta) fuse_tidy <-tidy(fuse_beta) edge2_tidy <-tidy(edge2_beta)# Get pseudo R² iucn_r2 <-cor(predict(iucn_beta), model_data$IUCN_protection_beta)^2 fuse_r2 <-cor(predict(fuse_beta), model_data$FUSE_protection_beta)^2 edge2_r2 <-cor(predict(edge2_beta), model_data$EDGE2_protection_beta)^2# Calculate marginal effects (average change in percentage per unit increase) iucn_margins <-margins(iucn_beta) fuse_margins <-margins(fuse_beta) edge2_margins <-margins(edge2_beta)# Extract AME and SE iucn_ame_summary <-summary(iucn_margins) fuse_ame_summary <-summary(fuse_margins) edge2_ame_summary <-summary(edge2_margins)# Convert to percentage (margins gives proportions) iucn_ame <- iucn_ame_summary$AME *100 iucn_ame_se <- iucn_ame_summary$SE *100 fuse_ame <- fuse_ame_summary$AME *100 fuse_ame_se <- fuse_ame_summary$SE *100 edge2_ame <- edge2_ame_summary$AME *100 edge2_ame_se <- edge2_ame_summary$SE *100cat("\n=== Average Marginal Effects (as percentages) ===\n")cat(paste("IUCN AME:", round(iucn_ame, 2), "%, SE:", round(iucn_ame_se, 2), "%\n"))cat(paste("FUSE AME:", round(fuse_ame, 2), "%, SE:", round(fuse_ame_se, 2), "%\n"))cat(paste("EDGE2 AME:", round(edge2_ame, 2), "%, SE:", round(edge2_ame_se, 2), "%\n"))# Calculate predicted values and differences# IUCN predictions iucn_pred <-predict(iucn_beta,newdata =data.frame(IUCN_status =1:5),type ="response") *100cat("\nPredicted protection by IUCN category:\n")print(data.frame(Category =1:5,IUCN_Category =c("LC", "NT", "VU", "EN", "CR"),Predicted_Protection =round(iucn_pred, 2))) iucn_avg_increase <-mean(diff(iucn_pred))cat(paste("Average increase per IUCN category:", round(iucn_avg_increase, 2), "%\n"))# FUSE predictions fuse_pred <-predict(fuse_beta,newdata =data.frame(FUSE_score =1:5),type ="response") *100cat("\nPredicted protection by FUSE score:\n")print(data.frame(Score =1:5,Predicted_Protection =round(fuse_pred, 2))) fuse_avg_increase <-mean(diff(fuse_pred))cat(paste("Average increase per FUSE unit:", round(fuse_avg_increase, 2), "%\n"))# EDGE2 predictions edge2_pred <-predict(edge2_beta,newdata =data.frame(EDGE2_score =1:5),type ="response") *100cat("\nPredicted protection by EDGE2 score:\n")print(data.frame(Score =1:5,Predicted_Protection =round(edge2_pred, 2))) edge2_avg_increase <-mean(diff(edge2_pred))cat(paste("Average increase per EDGE2 unit:", round(edge2_avg_increase, 2), "%\n"))# Create comprehensive summary table with AME comprehensive_summary <-data.frame(Approach =c("IUCN", "FUSE", "EDGE2"),Coefficient_logit =c( iucn_tidy$estimate[2], fuse_tidy$estimate[2], edge2_tidy$estimate[2] ),SE =c( iucn_tidy$std.error[2], fuse_tidy$std.error[2], edge2_tidy$std.error[2] ),z_value =c( iucn_tidy$statistic[2], fuse_tidy$statistic[2], edge2_tidy$statistic[2] ),P_value =c( iucn_tidy$p.value[2], fuse_tidy$p.value[2], edge2_tidy$p.value[2] ),AME_percent =c( iucn_ame, fuse_ame, edge2_ame ),AME_SE =c( iucn_ame_se, fuse_ame_se, edge2_ame_se ),Pseudo_R2 =c(iucn_r2, fuse_r2, edge2_r2) )cat("\n=== Comprehensive Beta Regression Summary ===\n")print(comprehensive_summary)# Save resultswrite.csv(comprehensive_summary, here::here("Outputs", "beta_regression_summary.csv"),row.names =FALSE)} else {cat("ERROR: Datasets have different number of rows.\n")cat("Protected fractions:", nrow(protected_fractions), "rows\n")cat("Shark metrics:", nrow(shark_metrics), "rows\n")}
Successfully joined datasets with 1000 species
=== IUCN Beta Regression ===
Call:
betareg(formula = IUCN_protection_beta ~ IUCN_status, data = model_data)
Quantile residuals:
Min 1Q Median 3Q Max
-2.5313 -0.5004 -0.0020 0.5608 2.0202
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.05978 0.07511 -0.796 0.426
IUCN_status 0.11768 0.02988 3.938 8.21e-05 ***
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 0.91473 0.03196 28.63 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 379.5 on 3 Df
Pseudo R-squared: 0.01504
Number of iterations: 7 (BFGS) + 1 (Fisher scoring)
=== FUSE Beta Regression ===
Call:
betareg(formula = FUSE_protection_beta ~ FUSE_score, data = model_data)
Quantile residuals:
Min 1Q Median 3Q Max
-2.3375 -0.4546 0.0051 0.5057 2.3049
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.19664 0.08966 -2.193 0.0283 *
FUSE_score 0.14261 0.06830 2.088 0.0368 *
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 1.04726 0.03673 28.52 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 211.9 on 3 Df
Pseudo R-squared: 0.004321
Number of iterations: 7 (BFGS) + 2 (Fisher scoring)
=== EDGE2 Beta Regression ===
Call:
betareg(formula = EDGE2_protection_beta ~ EDGE2_score, data = model_data)
Quantile residuals:
Min 1Q Median 3Q Max
-2.5490 -0.4954 -0.0052 0.4548 2.1662
Coefficients (mean model with logit link):
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.68964 0.09615 -7.173 7.34e-13 ***
EDGE2_score 0.46406 0.07688 6.036 1.58e-09 ***
Phi coefficients (precision model with identity link):
Estimate Std. Error z value Pr(>|z|)
(phi) 0.84174 0.02926 28.77 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Type of estimator: ML (maximum likelihood)
Log-likelihood: 506.7 on 3 Df
Pseudo R-squared: 0.04761
Number of iterations: 13 (BFGS) + 1 (Fisher scoring)
=== Average Marginal Effects (as percentages) ===
IUCN AME: 2.9 %, SE: 0.73 %
FUSE AME: 3.56 %, SE: 1.7 %
EDGE2 AME: 11.39 %, SE: 1.85 %
Predicted protection by IUCN category:
Category IUCN_Category Predicted_Protection
1 1 LC 51.45
2 2 NT 54.38
3 3 VU 57.28
4 4 EN 60.13
5 5 CR 62.92
Average increase per IUCN category: 2.87 %
Predicted protection by FUSE score:
Score Predicted_Protection
1 1 48.65
2 2 52.21
3 3 55.75
4 4 59.24
5 5 62.63
Average increase per FUSE unit: 3.5 %
Predicted protection by EDGE2 score:
Score Predicted_Protection
1 1 44.38
2 2 55.93
3 3 66.88
4 4 76.25
5 5 83.63
Average increase per EDGE2 unit: 9.81 %
=== Comprehensive Beta Regression Summary ===
Approach Coefficient_logit SE z_value P_value
IUCN_status IUCN 0.1176756 0.02987957 3.938332 8.205012e-05
FUSE_score FUSE 0.1426108 0.06830262 2.087925 3.680458e-02
EDGE2_score EDGE2 0.4640609 0.07687772 6.036350 1.576390e-09
AME_percent AME_SE Pseudo_R2
IUCN_status 2.898230 0.7271499 0.01750165
FUSE_score 3.558461 1.6986270 0.00239858
EDGE2_score 11.386873 1.8514148 0.02982912
Code
# Display the saved image in the Quarto documentknitr::include_graphics(here::here("Outputs", "Fig.2_combined_protection_dotplots.png"))
Manuscript maps
Individual maps
Code
CAPTAIN2_EDGE2_msmap <-ggplot() +geom_sf(data = CAPTAIN2_EDGE2_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill =NA, color ="lightgrey", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority EDGE2",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Global Conservation Priorities",#subtitle = "CAPTAIN2 - EDGE2 Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_EDGE2CAPTAIN2_FUSE_msmap <-ggplot() +geom_sf(data = CAPTAIN2_FUSE_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_FUSE, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_FUSE, fill =NA, color ="lightgrey", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority FUSE",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Global Conservation Priorities",#subtitle = "CAPTAIN2 - FUSE Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_FUSECAPTAIN2_IUCN_msmap <-ggplot() +geom_sf(data = CAPTAIN2_IUCN_sf, aes(color = Priority), size =0.5, alpha =0.7) +geom_sf(data = world_projected_CAPTAIN2_IUCN, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_IUCN, fill =NA, color ="lightgrey", size =0.5) +scale_color_gradientn(colors =c("white", "yellow", "darkblue"),values =c(0, 0.5, 1),name ="Priority IUCN",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Global Conservation Priorities",#subtitle = "CAPTAIN2 - IUCN Index, Budget: 0.1, Replicates: 50",x =NULL, y =NULL) + my_theme_CAPTAIN2_IUCN# Create a function to add labels (A), (B), etc.add_panel_labels <-function(plot, label) { plot +theme(plot.title =element_text(face ="bold", hjust =0, size =12) ) +labs(title =paste0("(", label, ")"))}# Add labels to each plot# First gridCAPTAIN2_IUCN_msmap_labeled <-add_panel_labels(CAPTAIN2_IUCN_msmap, "A")CAPTAIN2_FUSE_msmap_labeled <-add_panel_labels(CAPTAIN2_FUSE_msmap, "B")CAPTAIN2_EDGE2_msmap_labeled <-add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")# Combine plots into two separate grids, each with 3 rows and 1 columngrid1 <- CAPTAIN2_IUCN_msmap_labeled / CAPTAIN2_FUSE_msmap_labeled / CAPTAIN2_EDGE2_msmap_labeled# Save the grids if neededggsave(here::here("Outputs", "grid1_maps_continental_2ndrun_new.png"), grid1, width =8, height =15, dpi =300)
Difference maps
Code
# 1. IUCN - FUSE Difference MapIUCN_FUSE_msmap <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_FUSE_diff, aes(color = IUCN_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="lightgrey", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority (IUCN - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Difference in Conservation Priorities",#subtitle = "IUCN Index minus FUSE Index",x =NULL, y =NULL) + my_theme# 2. IUCN - EDGE2 Difference MapIUCN_EDGE2_msmap <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = IUCN_EDGE2_diff, aes(color = IUCN_minus_EDGE2), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="lightgrey", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority (IUCN - EDGE2)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Difference in Conservation Priorities",#subtitle = "IUCN Index minus EDGE2 Index",x =NULL, y =NULL) + my_theme# 3. EDGE2 - FUSE Difference MapEDGE2_FUSE_msmap <-ggplot() +geom_sf(data = globe_border, fill ="#F8F8F8", color =NA) +geom_sf(data = EDGE2_FUSE_diff, aes(color = EDGE2_minus_FUSE), size =0.5) +geom_sf(data = world_projected, fill ="lightgrey", color ="darkgrey", size =0.1) +geom_sf(data = globe_border, fill =NA, color ="lightgrey", size =0.5) +scale_color_gradientn(colors = diff_colors,limits =c(-1, 1),breaks =seq(-1, 1, by =0.25),labels =as.character(seq(-1, 1, by =0.25)),name ="Difference in Priority (EDGE2 - FUSE)",guide =guide_colorbar(barwidth =20, barheight =0.5, title.position ="top", title.hjust =0.5) ) +labs(#title = "Difference in Conservation Priorities",#subtitle = "EDGE2 Index minus FUSE Index",x =NULL, y =NULL) + my_theme# Create a function to add labels (A), (B), etc.add_panel_labels <-function(plot, label) { plot +theme(plot.title =element_text(face ="bold", hjust =0, size =12) ) +labs(title =paste0("(", label, ")"))}# Add labels to each plot# First gridIUCN_FUSE_msmap_labeled <-add_panel_labels(IUCN_FUSE_msmap, "A")IUCN_EDGE2_msmap_labeled <-add_panel_labels(IUCN_EDGE2_msmap, "B")EDGE2_FUSE_msmap_labeled <-add_panel_labels(EDGE2_FUSE_msmap, "C")# Combine plots into two separate grids, each with 3 rows and 1 columngrid2 <- IUCN_FUSE_msmap_labeled / IUCN_EDGE2_msmap_labeled / EDGE2_FUSE_msmap_labeled# Save the grids if neededggsave(here::here("Outputs", "grid3_maps_continental_2ndrun.png"), grid2, width =16, height =15, dpi =300)
Combined individual and difference maps
Code
library(patchwork)# Add labels to each plot# First gridCAPTAIN2_IUCN_msmap_labeled <-add_panel_labels(CAPTAIN2_IUCN_msmap, "A")CAPTAIN2_FUSE_msmap_labeled <-add_panel_labels(CAPTAIN2_FUSE_msmap, "B")CAPTAIN2_EDGE2_msmap_labeled <-add_panel_labels(CAPTAIN2_EDGE2_msmap, "C")# Second gridIUCN_FUSE_msmap_labeled <-add_panel_labels(IUCN_FUSE_msmap, "D")IUCN_EDGE2_msmap_labeled <-add_panel_labels(IUCN_EDGE2_msmap, "E")EDGE2_FUSE_msmap_labeled <-add_panel_labels(EDGE2_FUSE_msmap, "F")# Create each grid (3 rows, 1 column)grid1 <- CAPTAIN2_IUCN_msmap_labeled / CAPTAIN2_FUSE_msmap_labeled / CAPTAIN2_EDGE2_msmap_labeledgrid2 <- IUCN_FUSE_msmap_labeled / IUCN_EDGE2_msmap_labeled / EDGE2_FUSE_msmap_labeled# Combine the two grids side by side (2 columns)combined_grid <- grid1 | grid2# If you want to save itggsave(filename = here::here("Outputs", "Fig.3_combined_priority_maps.png"),plot = combined_grid,width =10, # Adjust width as needed for two columnsheight =12, # Adjust height as needed for three rowsdpi =300,bg ="white")# Display the saved image in the Quarto documentknitr::include_graphics(here::here("Outputs", "Fig.3_combined_priority_maps.png"))
Spatial distribution analysis
Code
# Filter for high priority cells (>0.9) and calculate percentages# Calculate total coastal cells and high-priority cells for each index# IUCNCAPTAIN2_IUCN_high <- CAPTAIN2_IUCN_sf %>%filter(Priority >0.9)iucn_total_cells <-nrow(CAPTAIN2_IUCN_sf)iucn_high_cells <-nrow(CAPTAIN2_IUCN_high)iucn_percent <-round((iucn_high_cells / iucn_total_cells) *100, 1)# FUSECAPTAIN2_FUSE_high <- CAPTAIN2_FUSE_sf %>%filter(Priority >0.9)fuse_total_cells <-nrow(CAPTAIN2_FUSE_sf)fuse_high_cells <-nrow(CAPTAIN2_FUSE_high)fuse_percent <-round((fuse_high_cells / fuse_total_cells) *100, 1)# EDGE2CAPTAIN2_EDGE2_high <- CAPTAIN2_EDGE2_sf %>%filter(Priority >0.9)edge2_total_cells <-nrow(CAPTAIN2_EDGE2_sf)edge2_high_cells <-nrow(CAPTAIN2_EDGE2_high)edge2_percent <-round((edge2_high_cells / edge2_total_cells) *100, 1)#Mean nearest neighbor for high priority grid cells:# Function to calculate Mean Nearest Neighbor Distancecalculate_mnn <-function(sf_high_priority) {# Get centroids and convert to coordinates coords <-st_coordinates(st_centroid(sf_high_priority))# Create a bounding box for the window bbox <-st_bbox(sf_high_priority)# Create observation window window <-owin(xrange =c(bbox["xmin"], bbox["xmax"]),yrange =c(bbox["ymin"], bbox["ymax"]))# Create point pattern object ppp_obj <-ppp(coords[,1], coords[,2], window = window)# Calculate mean nearest neighbor distance mnn <-mean(nndist(ppp_obj))# Calculate standard deviation for context mnn_sd <-sd(nndist(ppp_obj))return(list(mean = mnn, sd = mnn_sd))}# Calculate for each indexiucn_mnn <-calculate_mnn(CAPTAIN2_IUCN_high)fuse_mnn <-calculate_mnn(CAPTAIN2_FUSE_high)edge2_mnn <-calculate_mnn(CAPTAIN2_EDGE2_high)# Create summary data framemnn_summary <-data.frame(Index =c("IUCN", "FUSE", "EDGE2"),N_cells =c(iucn_high_cells, fuse_high_cells, edge2_high_cells),Percent =c(iucn_percent, fuse_percent, edge2_percent),MNN_mean =c(iucn_mnn$mean, fuse_mnn$mean, edge2_mnn$mean),MNN_sd =c(iucn_mnn$sd, fuse_mnn$sd, edge2_mnn$sd))print(mnn_summary)
# Method 1: If you have separate dataframes for each index# Assuming your data has columns: PUID, Priority, and geometry# Fxirst, identify high priority areas (>0.9) for each indexhigh_priority_EDGE2 <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$Priority >0.9, ]high_priority_FUSE <- CAPTAIN2_FUSE_sf[CAPTAIN2_FUSE_sf$Priority >0.9, ]high_priority_IUCN <- CAPTAIN2_IUCN_sf[CAPTAIN2_IUCN_sf$Priority >0.9, ]# Find congruent areas (present in all three)# Method using PUID (assuming you have PUID column)congruent_PUIDs <-intersect(intersect(high_priority_EDGE2$PUID, high_priority_FUSE$PUID), high_priority_IUCN$PUID)# Extract congruent areascongruent_areas <- CAPTAIN2_EDGE2_sf[CAPTAIN2_EDGE2_sf$PUID %in% congruent_PUIDs, ]# Method 2: If you need to merge dataframes first# Create a combined dataframe with all three priority scores# First, extract just the data (without geometry) from the other sf objectsFUSE_data <-st_drop_geometry(CAPTAIN2_FUSE_sf[, c("PUID", "Priority")])IUCN_data <-st_drop_geometry(CAPTAIN2_IUCN_sf[, c("PUID", "Priority")])# Merge with the EDGE2 sf object (keeping geometry)combined_priorities <-merge(CAPTAIN2_EDGE2_sf, FUSE_data, by ="PUID", suffixes =c("_EDGE2", "_FUSE"))combined_priorities <-merge(combined_priorities, IUCN_data, by ="PUID")names(combined_priorities)[names(combined_priorities) =="Priority"] <-"Priority_IUCN"# Identify congruent high priority areascongruent_areas_v2 <- combined_priorities[combined_priorities$Priority_EDGE2 >0.9& combined_priorities$Priority_FUSE >0.9& combined_priorities$Priority_IUCN >0.9, ]# Create a map showing only congruent areascongruent_map <-ggplot() +geom_sf(data = congruent_areas, aes(fill ="Congruent High Priority"), color ="red", size =0.8, alpha =0.8) +geom_sf(data = world_projected_CAPTAIN2_EDGE2, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_CAPTAIN2_EDGE2, fill =NA, color ="black", size =0.5) +scale_fill_manual(values =c("Congruent high priority"="red"),name ="Priority areas") +labs(title ="Congruent high priority areas",subtitle ="Areas with priority > 0.9 in IUCN, FUSE and EDGE2 dimensions",x =NULL, y =NULL) + my_theme_CAPTAIN2_EDGE2# Display the mapprint(congruent_map)
Code
ggsave(here::here("Outputs", "congruent_high_prioritiy_areas_09.png"), congruent_map, width =10, height =8, dpi =300, bg ="white")# Create detailed congruence categories# Define logical conditions for each indexEDGE2_high <- combined_priorities$Priority_EDGE2 >0.9FUSE_high <- combined_priorities$Priority_FUSE >0.9IUCN_high <- combined_priorities$Priority_IUCN >0.9# Create detailed congruence categoriescombined_priorities$congruence_category <-case_when( EDGE2_high & FUSE_high & IUCN_high ~"All three indices",!EDGE2_high & FUSE_high & IUCN_high ~"IUCN + FUSE", # Changed from "FUSE + IUCN" EDGE2_high &!FUSE_high & IUCN_high ~"IUCN + EDGE2", # Changed from "EDGE2 + IUCN" EDGE2_high & FUSE_high &!IUCN_high ~"FUSE + EDGE2", # Changed from "EDGE2 + FUSE"!EDGE2_high &!FUSE_high & IUCN_high ~"IUCN only",!EDGE2_high & FUSE_high &!IUCN_high ~"FUSE only", EDGE2_high &!FUSE_high &!IUCN_high ~"EDGE2 only",TRUE~"No high priority")# Convert to factor with desired ordercombined_priorities$congruence_category <-factor( combined_priorities$congruence_category,levels =c("All three indices","IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2", # Changed order"IUCN only", "FUSE only", "EDGE2 only"))# Filter data to only include high priority areashigh_priority_data <- combined_priorities[combined_priorities$congruence_category %in%c("All three indices", "IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2", # Changed names"IUCN only", "FUSE only", "EDGE2 only"), ]# If simple plot fails, let's check the data more thoroughlyif (nrow(high_priority_data) ==0) {cat("No high priority data found! Checking original data...\n")cat("EDGE2 > 0.9:", sum(combined_priorities$Priority_EDGE2 >0.9, na.rm =TRUE), "\n")cat("FUSE > 0.9:", sum(combined_priorities$Priority_FUSE >0.9, na.rm =TRUE), "\n") cat("IUCN > 0.9:", sum(combined_priorities$Priority_IUCN >0.9, na.rm =TRUE), "\n")}# Extract coordinates for the congruence plotif (nrow(high_priority_data) >0) { coords <-st_coordinates(st_centroid(high_priority_data)) plot_data <-data.frame(x = coords[,1],y = coords[,2], category = high_priority_data$congruence_category )# Remove any rows with NA category plot_data <- plot_data[!is.na(plot_data$category), ]#cat("Creating enhanced congruence plot...\n")#cat("Plot data dimensions:", nrow(plot_data), "points\n")# Create the enhanced congruence plot congruence_map <-ggplot(plot_data, aes(x = x, y = y, color = category)) +geom_point(size =1.2, alpha =0.85, stroke =0) +# Larger points, no stroke to reduce overlapscale_color_manual(values =c("All three indices"="#8B0000", # Dark red"EDGE2 + FUSE"="#FF8C00", # Dark orange "EDGE2 + IUCN"="#DC143C", # Crimson"FUSE + IUCN"="#FFD700", # Gold"EDGE2 only"="#4169E1", # Royal blue"FUSE only"="#32CD32", # Lime green"IUCN only"="#9370DB"), # Medium purplename ="High Priority\nCongruence",guide =guide_legend(override.aes =list(size =4, alpha =1), # Larger, more opaque legend pointstitle.position ="top",title.hjust =0.5,ncol =4# Two columns for legend to save space )) +labs(title ="Global Conservation Priority Congruence Analysis",subtitle ="High priority areas (>0.9) showing agreement patterns across EDGE2, FUSE, and IUCN indices",x ="Longitude", y ="Latitude",caption =paste("Total high priority areas:", nrow(plot_data)) ) +theme_minimal() +theme(# Plot aestheticsplot.background =element_rect(fill ="white", color =NA),panel.background =element_rect(fill ="#f8f9fa", color =NA),panel.grid.major =element_line(color ="white", size =0.5, linetype ="solid"),panel.grid.minor =element_line(color ="white", size =0.25, linetype ="solid"),# Text stylingplot.title =element_text(size =16, hjust =0.5, face ="bold", margin =margin(b =10)),plot.subtitle =element_text(size =12, hjust =0.5, color ="gray30",margin =margin(b =20)),plot.caption =element_text(size =10, color ="gray50", hjust =1),# Axis stylingaxis.title =element_text(size =12, face ="bold"),axis.text =element_text(size =10, color ="gray30"),axis.ticks =element_line(color ="gray50", size =0.5),# Legend stylinglegend.position ="bottom",legend.background =element_rect(fill ="white", color ="gray80", size =0.5),legend.margin =margin(t =15, r =10, b =10, l =10),legend.title =element_text(size =12, face ="bold"),legend.text =element_text(size =10),legend.key.size =unit(0.8, "cm"),# Panel borderpanel.border =element_rect(color ="gray70", fill =NA, size =0.5) ) +# Add coordinate system and aspect ratiocoord_fixed(ratio =1.3) +# Adjust ratio for better map appearance# Add subtle borders around plot areascale_x_continuous(expand =c(0.02, 0.02)) +scale_y_continuous(expand =c(0.02, 0.02))# Print the enhanced plot#print(congruence_map)#cat("Enhanced congruence plot created successfully!\n")# Print summary by categorycat("\nBreakdown by congruence category:\n") category_counts <-table(plot_data$category) category_percentages <-round(prop.table(category_counts) *100, 1)for(i in1:length(category_counts)) {cat(sprintf("%-20s: %4d points (%4.2f%%)\n", names(category_counts)[i], category_counts[i], category_percentages[i])) }} else {cat("No valid high priority data to plot\n")}
Breakdown by congruence category:
All three indices : 1459 points (46.50%)
IUCN + FUSE : 75 points (2.40%)
IUCN + EDGE2 : 1130 points (36.00%)
FUSE + EDGE2 : 275 points (8.80%)
IUCN only : 29 points (0.90%)
FUSE only : 17 points (0.50%)
EDGE2 only : 154 points (4.90%)
Code
# Print summary of congruence patterns#cat("\nCongruence Pattern Summary:\n")congruence_summary <-table(combined_priorities$congruence_category, useNA ="ifany")#print(congruence_summary)# Calculate percentagestotal_high_priority <-sum(congruence_summary[names(congruence_summary) !="No high priority"], na.rm =TRUE)congruence_percentages <-round(congruence_summary / total_high_priority *100, 2)# ms style map -----# Step 1: Check original data#print("=== STEP 1: CHECKING ORIGINAL DATA ===")#print(head(plot_data))#print(summary(plot_data[c("x", "y")]))#print(unique(plot_data$category))# Step 2: Transform to sf objectcongruence_sf <-st_as_sf( plot_data, coords =c("x", "y"), crs = mcbryde_thomas_2 # Assuming WGS84)# Step 3: Project to McBryde-Thomas 2congruence_sf <-st_transform(congruence_sf, crs = mcbryde_thomas_2)# Step 4: Check world projectionworld_projected_congruence <-st_transform(world, crs = mcbryde_thomas_2)# Step 5: Create globe borderglobe_bbox <-rbind(c(-180, -90), c(-180, 90), c(180, 90), c(180, -90), c(-180, -90))globe_border_congruence <-st_polygon(list(globe_bbox)) %>%st_sfc(crs =4326) %>%st_sf(data.frame(rgn ='globe', geom = .)) %>% smoothr::densify(max_distance =0.5) %>%st_transform(crs = mcbryde_thomas_2)# Step 6: Check intersectionpoints_in_globe <-st_intersects(congruence_sf, globe_border_congruence, sparse =FALSE)# Step 7: Extract coordinates for verificationcoords <-st_coordinates(congruence_sf)# Step 8: Create thememy_theme_congruence <-theme_minimal() +theme(legend.position ="bottom",legend.direction ="horizontal",legend.box ="vertical",legend.margin =margin(t =20, r =0, b =0, l =0),legend.title =element_text(margin =margin(b =10)),panel.background =element_rect(fill ="white", color =NA),plot.background =element_rect(fill ="white", color =NA),panel.grid =element_blank(),axis.text =element_blank(), # Add this lineaxis.ticks =element_blank() # Add this line )# Step 9: Create plot with debuggingcongruence_map <-ggplot() +geom_sf(data = congruence_sf, aes(color = category), size =1.2, alpha =0.85) +geom_sf(data = world_projected_congruence, fill ="lightgrey", color ="lightgrey", size =0.1) +geom_sf(data = globe_border_congruence, fill =NA, color ="lightgrey", size =0.5) +scale_color_manual(values =c("All three indices"="#8B0000","IUCN + FUSE"="#FFD700", # Changed from "FUSE + IUCN""IUCN + EDGE2"="#DC143C", # Changed from "EDGE2 + IUCN""FUSE + EDGE2"="#FF8C00", # Changed from "EDGE2 + FUSE""IUCN only"="#9370DB","FUSE only"="#32CD32","EDGE2 only"="#4169E1"),breaks =c("All three indices","IUCN + FUSE", "IUCN + EDGE2", "FUSE + EDGE2","IUCN only", "FUSE only", "EDGE2 only"), # Finally: FUSE + EDGE2name ="High priority congruence",guide =guide_legend(override.aes =list(size =4, alpha =1),title.position ="top", title.hjust =0.5,ncol =4 )) +labs(title ="Global conservation priority congruence analysis",subtitle ="High priority areas (>0.9) showing agreement patterns across IUCN, FUSE and EDGE2 prioritisations",x =NULL, y =NULL ) + my_theme_congruenceggsave(here::here("outputs", "congruent_high_prioritiy_areas_09_all_indices.png"), congruence_map, width =10, height =8, dpi =300, bg ="white")# Display both plotsprint(congruence_map)